Porous Media Flow

Flow through porous media is governed by Darcy’s law

Source left: https://serc.carleton.edu/integrate/teaching_materials/water_science_society/student_materials/926 Source right: https://pressbooks.cuny.edu/gorokhovich/chapter/porosity-and-permeability-darcy-law/

Needed: Technique to determine macroscale process

Source left: Cross-section through a fiber-reinforced composite with two orientations

(University of Augsburg)

Let’s start simple: Solve directly

Stationary diffusion at fixed boundary conditions:

\[ \begin{aligned} \frac{d}{dx} \left( a^\epsilon(x) \frac{d}{dx} u^\epsilon (x) \right) & = 0 \qquad \text{for} \quad 0 \leq x \leq L \end{aligned} \]

with an oscillating model parameter

\[ \begin{aligned} a^\epsilon(x) &= \frac{1}{1 + 2 \sin^2 (\pi x / \epsilon)} \end{aligned} \]

Consider the properties of the parameter

We are interested in the limit \(\epsilon \Rightarrow 0\)!

Challenge: The parameter does not converge (pointwise)!

Solution: Solve for a sequence of parameters

Observation: The numerical solution seems to converge!

We postulate \(\lim_{\epsilon \rightarrow 0} u^\epsilon (x) = x\)!

Mathematical homogenization

Idea 1: Identify macro-scale \(L\), micro-scale \(l\) and introduce scaling

\[ \begin{equation*} y := \frac{L}{l} x = \frac{x}{\epsilon} \qquad \text{in which} \quad \epsilon := \frac{l}{L} \end{equation*} \]

Idea 2: Interpret functions as two-scale and make asymptotic expansion Ansatz:

\[ \begin{equation*} f^\epsilon = f(x,y) = \sum_{i=0}^n \epsilon^i f^{(i)}(x,y) = f^{(0)}(x,y) + \epsilon^1 f^{(1)}(x,y) + \epsilon^2 f^{(2)}(x,y) + \mathcal O (\epsilon^3) \end{equation*} \]

Differentiation with respect to \(x\) then translates according to:

\[ \begin{aligned} \frac{d}{dx} u(x,y) &= \frac{d}{dx} u\left(x,y(x)\right) \\ &= \partial_x u(x,y) + \partial_y u(x,y) \underbrace{\frac{d}{dx}y(x)}_{= \epsilon^{-1}} = \partial_x u + \epsilon^{-1} \partial_y u \end{aligned} \]

Mathematical homogenization

Substitution yields

\[ \begin{align*} & \quad \frac{d}{dx} \left( a(y) \frac{d}{dx} \left(u^{(0)} + \epsilon u^{(1)} + \epsilon^2 u^{(2)}\right) \right) \\ & = \frac{d}{dx} \left( a(y) \left(\partial_x u^{(0)} + \epsilon^{-1} \partial_y u^{(0)} + \epsilon \partial_x u^{(1)} + \partial_y u^{(1)} + \epsilon^2 \partial_x u^{(2)} + \epsilon \partial_y u^{(2)}\right) \right)\\ % first line & = \underbrace{\partial_x \left( a(y) \, \partial_x u^{(0)} \right)}_{\text{order} \; \epsilon^0} + \underbrace{\epsilon^{-1}\,\partial_y \left( a(y) \, \partial_x u^{(0)} \right)}_{\text{order} \; \epsilon^{-1}} + \underbrace{\epsilon^{-1}\,\partial_x \left( a(y) \, \partial_y u^{(0)} \right)}_{\text{order} \; \epsilon^{-1}} + \underbrace{\epsilon^{-2}\,\partial_y \left( a(y) \, \partial_y u^{(0)} \right)}_{\text{order} \; \epsilon^{-2}} +\\ % second line & + \underbrace{\epsilon \, \partial_x \left( a(y) \, \partial_x u^{(1)} \right)}_{\text{order} \; \epsilon^1} + \underbrace{\partial_y \left( a(y) \, \partial_x u^{(1)} \right)}_{\text{order} \; \epsilon^{0}} + \underbrace{\partial_x \left( a(y) \, \partial_y u^{(1)} \right)}_{\text{order} \; \epsilon^{0}} + \underbrace{\epsilon^{-1}\,\partial_y \left( a(y) \, \partial_y u^{(1)} \right)}_{\text{order} \; \epsilon^{-1}} +\\ % third line & + \underbrace{\epsilon^2 \, \partial_x \left( a(y) \, \partial_x u^{(2)} \right)}_{\text{order} \; \epsilon^2} + \underbrace{\epsilon \,\partial_y \left( a(y) \, \partial_x u^{(2)} \right)}_{\text{order} \; \epsilon^{1}} + \underbrace{\epsilon \,\partial_x \left( a(y) \, \partial_y u^{(2)} \right)}_{\text{order} \; \epsilon^{1}} + \underbrace{\partial_y \left( a(y) \, \partial_y u^{(2)} \right)}_{\text{order} \; \epsilon^{0}} + \end{align*} \]

Mathematical homogenization

Re-ordering in powers of \(\epsilon\) yields

\[ \begin{align*} \mathcal O \left(\epsilon^{-2}\right): \quad &\partial_y \left( a(y) \, \partial_y u^{(0)} \right) = 0 \\ \mathcal O \left(\epsilon^{-1}\right): \quad &\partial_x \left( a(y) \, \partial_y u^{(0)} \right) + \partial_y \left( a(y) \, \partial_x u^{(0)} \right) + \partial_y \left( a(y) \, \partial_y u^{(1)} \right) = 0\\ \mathcal O \left(\epsilon^{0}\right): \quad &\partial_x \left( a(y) \, \partial_x u^{(0)} \right) + \partial_y \left( a(y) \, \partial_x u^{(1)} \right) + \partial_x \left( a(y) \, \partial_y u^{(1)} \right) + \partial_y \left( a(y) \, \partial_y u^{(2)} \right) = 0\\ \mathcal O \left(\epsilon^{1}\right): \quad &\partial_x \left( a(y) \, \partial_x u^{(1)} \right) + \partial_y \left( a(y) \, \partial_x u^{(2)} \right) + \partial_x \left( a(y) \, \partial_y u^{(2)} \right) + \partial_y \left( a(y) \, \partial_y u^{(3)} \right) = 0 \end{align*} \]

Boundary conditions on both scales are

\[ \begin{align*} &u^{(0)} (0,y) = 0 \qquad u^{(0)} (1,y) = 1 \\ &u^{(i)} (0,y) = 0 \qquad u^{(i)} (1,y) = 0 \qquad &\text{for all} \; i \geq 1 \\ &u^{(i)} (x,0) = u^{(i)} (x,1) \qquad &\text{for all} \; i \geq 0 \\ &\partial_y u^{(i)} (0,y) = \partial_y u^{(i)} (1,y) = 0 \qquad &\text{for all} \; i \geq 0 \end{align*} \]

Solve step by step

Step 1: We observe that \(u^{(0)}\) independent of \(y\) is a solution to the largest scale:

\[ u^{(0)}(x,y) = u^{(0)}(x) \]

Step 2: Substituting the latter into the \(\epsilon^{-1}\) yields \[ \begin{align*} \underbrace{\partial_y \left( a(y) \, \partial_x u^{(0)} \right)}_{\left( \partial_y a(y) \right) \partial_x u^{(0)} + \underbrace{a(y) \partial_{yx} u^{(0)}}_{=0}} + \partial_y \left( a(y) \, \partial_y u^{(1)} \right) = 0. \end{align*} \]

Note, that with \(u^{(1)}\) defined according to \[ \begin{align*} u^{(1)}(x,y) = \left( \frac{\int_0^y \frac{1}{a(\tilde y)} d\tilde y}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} - y + c\right) \partial_x u^{(0)} + f(x) \end{align*} \]

this equation is fulfilled.

Solve step by step

Step 3: Integrate the \(\epsilon\) scale over a periodic cell \(Y\): \[ \begin{align*} & \int_0^1 \partial_x \left( a(y) \partial_x u^{(0)} \right) d y + \underbrace{\int_0^1 \partial_y \left( a(y) \, \partial_x u^{(1)} \right) d y}_{[a(y) \partial_x u^{(1)}]_0^1 = 0} \\ & + \int_0^1 \partial_x \left( a(y) \partial_y u^{(1)} \right) d y + \underbrace{\int_0^1 \partial_y \left( a(y) \partial_y u^{(2)} \right) d y}_{[a(y) \partial_y u^{(2)}]_0^1 = 0} \end{align*} \]

where we observe several terms to cancel out due to the micro-scale periodicity.The following two terms remain: \[ \begin{align*} \int_0^1 \partial_x \left( a(y) \partial_x u^{(0)} \right) + \partial_x \left( a(y) \partial_y u^{(1)} \right) d y \end{align*} \]

Solve step by step

Substituting in \(u^{(1)}\) as defined before finally yields

\[ \begin{align*} & \int_0^1 \partial_x \left( a(y) \partial_x u^{(0)} \right) + \partial_x \left( a(y) \left( \left( \frac{\frac{1}{a( y)}}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} - 1 \right) \partial_x u^{(0)} \right) \right) d y \\ = & \int_0^1 \partial_x \left( \frac{1}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} \partial_x u^{(0)} \right) d y = \partial_x \left(\underbrace{\int_0^1 \frac{1}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} d y }_{=:a^{(0)}}\right) \partial_x u^{(0)} \\ = & \frac{d}{dx} \left( a^{(0)} \frac{d}{dx} u^{(0)} \right). \end{align*} \]

Solution

For our particular example, we get \(\int_0^1 \frac{1}{\int_0^1 \frac{1}{a(\tilde y)} d\tilde y} d y = \frac{1}{2}\)

and hence

\[ \begin{align*} u^{(0)} = x \qquad \text{and} \qquad u^{(1)} = - \frac{1}{ 4 \pi} \sin \left( 2 \pi y \right). \end{align*} \]

The complete solution reads \[ \begin{align*} u(x,y) = u^{(0)} + \epsilon u^{(1)} = x - \frac{\epsilon}{ 4 \pi} \sin \left( 2 \pi y \right) \end{align*} \]

and we indeed see that \[ \begin{align*} \underset{\epsilon \to 0}{\lim} \, u(x,y) = x, \end{align*} \]

as we initially postulated!

Homogenized solution

Homogenized solution corresponds to direct solution!

Postulated \(\lim_{\epsilon \rightarrow 0} u^\epsilon (x) = x\) has been proven!